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Abstract 

We obtain highly accurate solutions to the Thomas-Fermi equations for atoms and 
atoms in very strong magnetic fields. We apply the Pade-Hankel method, numerical 
integration, power series with Pade and Hermite-Pade approximants and Chebyshev 
polynomials. Both the slope at origin and the location of the right boundary in the 
magnetic-field case are given with unprecedented accuracy. 



1 Introduction 



The Thomas-Fermi model is one of the simplest approaches to the study of 
the potential and charge densities in a variety of systems, like, for example, 
atoms [HE], molecules [HE], atoms in strong magnetic fields [rj||gHTT]. met- 
als and crystals [121 [33] and dense plasmas [14]. For this reason there has 
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been great interest in the accurate calculation of the solution to the Thomas- 
Fermi equation [T|l3|l5| IT5H2~i] . In particular the accurate results obtained by 
Kobayashi et al [3] by numerical integration are commonly chosen as bench- 
mark data for testing other approaches. The even more accurate results of 
Rijnierse do not appear to be so well known, probably because they do not 
appear to have been published and are only quoted in the book by Torrens [5] . 

The behaviour of the solution to the nonlinear Thomas-Fermi equation de- 
pends on the slope at the origin. The critical slope at the origin suitable for 
neutral atoms is of particular interest and has been estimated by many au- 
thors (see, for example, Kobayashi et al [3]). Particularly accurate results for 
this critical slope were obtained by Amore and Fernandez [T9l[2D] and later 
by Fernandez [211 [22] using the Pade-Hankel method (PHM). Abbasbandy 
and Bervillier [23] considerably improved this estimate by means of a judi- 
cious analytic continuation of the expansion of the solution about the origin 
and Boyd [21] reported an even more accurate result obtained by means of a 
rational Chebyshev series. 

The analytic properties of the solution of the Thomas-Fermi equation under 
different boundary conditions (in addition to the physically relevant ones) 
are also of great interest and have been studied by several authors (see, for 
example, Hille [2HI2E] and the references therein). 

The purpose of this paper is twofold. First, we want to stress the different 
behaviour of the solution of the Thomas-Fermi equation for atoms and atoms 
in strong magnetic fields that was overlooked in a recent application of the 
PHM [22]- More precisely, in his application of the PHM Fernandez [22] as- 
sumed the incorrect asymptotic behaviour at infinity suggested by Banerjee et 
al [8]. However, since the PHM does not take into account the second (outer) 
boundary condition explicitly Fernandez obtained an accurate slope at the 
origin. Our second goal is to show that the available computer-algebra soft- 
ware enable one to solve the Thomas-Fermi equations (and, certainly, other 
nonlinear equations as well) with great accuracy. We want to provide suffi- 
ciently accurate solutions that may be used as benchmark data for testing 
future analytical or numerical methods. 

Boyd [21] obtained the most accurate critical slope at the origin of the solution 
to the Thomas-Fermi equation for neutral atoms by means of rational Cheby- 
shev series. In this paper we try the Chebyshev polynomials on the equation 
for neutral atoms in very strong magnetic fields. 

In Section [2] we outline the expansions of the solution to the Thomas- Fermi 
equation about origin, at infinity, about poles and zeroes/branch points. Such 
expansions have already been discussed by other authors and also used as aids 
for obtaining approximate analytical solutions as well as accurate numerical 
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results [IH3l|5H71[IMHl[23l|25l[271|28]. The main purpose of this section is to 
stress the difference between the Thomas-Fermi equations for isolated atoms 
and atoms in strong magnetic fields. In section [3] we describe the accurate 
results obtained by means of the PHM and from straightforward integration 
of the differential equations. In Sec. H]we discuss the application of two power 
series and their Pade and Hermite-Pade approximants to the Thomas-Fermi 
equation for a neutral atom in a magnetic field. In Sec. [5] we apply Chebyshev 
polynomials to the same proble and in Sec. 15.21 we consider analytical results 
based such polynomials of small degree. Finally, in Sec.Elwe draw conclusions. 



2 Expansions for the solutions to the Thomas-Fermi equations 



In order to facilitate the discussion and to make this paper clearer in this sec- 
tion we summarize the well known expansions of the Thomas- Fermi equations 
about some characteristic points. As indicated above, such expansions are 
well known and have been widely used by other authors for several different 
purposes [THSISHZIIIMHII2311251E3I28]. 

In this paper we restrict ourselves to the simplest cases. The first one is the 
Thomas-Fermi equation for an atom 



u{xf' 2 

«(0) = 1 (1) 

It can be expanded about origin as 

4x 3/2 2ax 5/2 x 3 , s 

u = l + ax + —— + — — + — + ... (2) 

3 5 3 



where a = u'(0) < is the unknown slope at the origin. 

There is a critical slope u' and the behaviour of the solution depends on the 
relation between u'(0) and u' Q . If u'(0) < u' the solution vanishes at a movable 
branch point [28J x = x q around which it behaves as 

u = ta + ^ /_ + _ + _+! (3) 



where s = (x q — x)/x q and b is a constant. 
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If u'(0) > u' the solution decreases and reaches a minimum at x = x m about 
which it behaves as 



u = a + (qqX 



3/2 

m ) 



1 h A h A '- h . . . (4) 

2 12 32 320 



where s = (x m — x)/x m and ao is a positive constant. 

After this minimum the solution increases and tends to infinity because of a 
movable singularity at x = x s about which it behaves as 

1 /400 2000 2000 10000 3189895 \ 
U ~x^\~7 4 9V" ~ 81^" ~ 729s ~ 177147 + " 7 ^ ' 



where s = (x s — x)/x s . Abbasbandy and Bervillier [23] derived an alternative 
and more detailed expansion in terms of the variable z = x 1//2 for the function 
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It is also well known that x q — > oo as u'(0) — > u' from the left and x m , x s — > oo 
as w'(0) — > u' Q from the right. At this limit the solution tends monotonically 
to zero according to 



u 



^ + dx(I -^ + !^+5! dV -^ + ... (6 ) 



where d is a constant. There have been reasonably successful attempts at 
matching the expansions (J2J and ([6]) by means of appropriate nonlinear trans- 
formations of the independent variable [T71IT8] . 

From a physical point of view the zero at x = x q (see the expansion (jSJ) is 
the radio of an atom if 

N 

- x q u(x q ) = l- — = q (7) 



where N is the number of electrons, Z is the atomic number and q is the degree 
of ionization (note that N < Z). For a neutral atom (q = 0) the boundary 
condition u(x ) = u\xq) = takes place at x = oo as indicated above for the 
case w'(0) = u' . 

Abbasbandy and Bervillier [23] argued that the PHM applies successfully to 
this problem because the Hankel condition sends the movable singularity at 
x = x s to infinity. 
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The Thomas- Fermi equation for an atom in a strong magnetic field is (Tomishina 
and Yonei [9] proposed a somewhat more realistic model that we do not discuss 
here) 



u"(x) = \Jxu(x) 

u(0) = l (8) 
In this case the expansion about the origin is given by 

Ax 5 / 2 2ax 7 ' 2 

u=l + ax + — — + — — + . . . 9 
15 35 

where a = u'(0). As in the preceding case the behaviour of the solution depends 
on this slope at the origin that also exhibits a critical value u' . 

If u'(0) < u' the solution vanishes at a movable branch point x = x q according 
to 



" = b ° + ± TE + < 10) 

where s = (x q — x)/x q and h is a constant. 

If u'(0) > u'q the solution exhibits a minimum at x = x m around which it 
behaves as 



v a X m 2 \J a X m 3 % X m ~~ \J a Q X m 4 ^>\J a X m + & X m 5 , , 

11 — CLn -\- s — s -\- s — s — I — . . . ( J L 

2 12 96 960 V ' 

where s = (x m — x)/x m and ao is a constant. 

The main difference between this equation and the preceding one is that in 
this case the pole is located at infinity. For large values of the coordinate the 
solution behaves as 

„ = ± + + 4 + f ) ^ + ■ ■ ■ (12) 



where Co is a constant. 

When u'(0) = u' Q the solution and its first derivative vanishes at x = xq 
according to the expansion 

5 4 ( 1 S S 2 gS \ 

u — xts ... (13) 

\ 144 336 7056 16464 ) v ; 
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where s = (xq — x)/xq. The location of the minimum x m approaches xq from 
below as u'(0) approaches u' from above. The zero/branch point x q also ap- 
proaches Xq from below as u'(0) approaches u' from below. However, the func- 
tion is analytic at x = Xo when tt'(O) = u' as shown in Eq. fll3p . The solution 
with the critical slope at the origin also tends to infinity as x — > oo according 
to equation ffl2|) . According to Banerjee et al [8] the universal solution cor- 
responding to neutral atoms (N = Z) satisfies u(xq) = and u'(xq) = at 
Xo = oo. Based on this conjecture Fernandez [22] applied the PHM and ob- 
tained the somewhat more accurate value of u' Q . However Hill et al [TT] showed 
that x is finite and can be related to x q by means of a perturbation expansion 
of the form. 



which clearly shows that x q approaches Xq from below as q — > (and u'(0) — > 
u' from below) . They confirmed this result by numerical integration of Eq. ([8]) . 

In the two cases discussed above the PHM yields the correct critical slope at 
the origin disregarding the second boundary condition. In the first example 
both u(x) and u'(x) vanish as x — > oo, while in the second example they vanish 
at a finite value Xq of the independent variable. Abbasbandy and Bervillier [25] 
suggested that the success of the PHM is based on "forcing the localization at 
infinity of a movable singularity (when it exists)". They also stated that "if the 
second boundary is located at infinity, the PHM has a particular significance" . 
The success of the PHM for the Thomas-Fermi equation shows that the 
approach is also suitable when the second boundary condition takes place at 
a finite point. In this case the movable singularity pushed to infinity may 
be the zero/branch point x q discussed above (Eq. ( fit)]) ). It was shown that x q 
approaches xo as u'(0) approaches the critical slope but it may jump to infinity 
when w'(0) = u' Q leaving the solution analytic at x as discussed above. We 
cannot prove this conjecture rigorously but we believe that it sounds plausible. 



3 PHM and numerical integration 

We first review the main points of the PHM. Following Amore and Fernandez 



[rgH22] we choose the new variables x = t 2 and v(t) = ^u{t 2 ). The function 

v (t) can be expanded in a Taylor series 



x q = x - (2Aq/x 2 ) 1/3 + ... 



(14) 




oo 



V(t)= 5>^ 



(15) 



j=0 



where the coefficients Vj, j > 4 depend on t> 2 = u' /2. 
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We construct the Hankel determinants Hfy = \vi + j + d+i\ i -_ 1 D _ x that depend 
on the unknown slope at the origin u' and obtain sequences of roots u' Q [D, d] 
of Hp = 0, D = 2, 3, . . ., (d = 0, 1, . . . fixed) that converge towards the critical 
slope. 

We carried out PHM calculations for values of D greater than those used before 
[22] and also numerical integration based commands built in Mathematica 
together with the bisection method. We describe the results in what follows. 

There are far too many results for the critical slope of the solution to the 
Thomas-Fermi equation for neutral atoms. Tabled] just shows the most accu- 
rate ones. Present PHM result was estimated from sequences of roots of the 
Hankel determinants with D < 60 and d = 3. 

The Thomas-Fermi equation for a neutral atom in a strong magnetic field 
has not been so widely studied. Table [2] shows the available critical slopes at 
origin. Present PHM result was estimated by comparing results with D < 60 
and d — 1, 3. The PHM exhibits a much greater rate of convergence for this 
problem. 

In order to provide benchmark data for testing other approaches in the future 
we have calculated u(x) and u'(x) for Eq. (OQ) and u(x) for Eq. ([5]) as accurately 
as possible using straightforward numerical integration. The results are shown 
in tables E] IH and 



4 Power-series approaches 

4-1 Power- series for the original equation 

When u'(0) = u' the function is analytic at and numerical experimen- 

tation suggests that the Taylor series f)13p converges over the whole domain 
of interest, even at the left endpoint x — 0. It is convenient to rescale the 
coordinate by defining s = (xq — x)/xq so that x = corresponds to s — — 1. 
Therefore we may try to calculate the parameters u' and Xq by means of the 
sequences of partial sums 

M X — X 

u M (x) = xls 4 ^2cjS j , s = (16) 

and the conditions ma/(0) = 1, u' M (0) = u' . Since the singularity at s = — 1 
(x = 0) is proportional to x 5//2 (see Eq. ([9])) the coefficients Cj decrease as 
0(j -7 / 2 ) and the error of the M — th partial sum falls proportionally to M~ 5//2 . 
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For example, in this way we obtain x$ ~ 3.06882 that is correct for five digits 
versus the actual value Xq ~ 3.068857. It is remarkable that a power series 
is able to yield any accuracy at all for a nonlinear, singular boundary value 
problem. 



4-2 Quadratic Pade approximants 



Accelerating the convergence by applying ordinary Pade approximations pro- 
duced no improvement because the function u has a branch point at the left 
boundary, precisely where we are summing the series to approximate u'(0). 
However, the so-called "Hermite-Pade" or "Shafer" approximation is much 
more successful. 

The quadratic Shafer approximant u[K/L/M\(s) is defined to be the solution 
of the quadratic equation [i8H5Tj 

P{s) {u[K/L/M}) 2 + Q{s) u[K/L/M] + R(s) = (17) 



where the polynomials P, Q and R are of degrees K, L and M, respec- 
tively. These polynomials are chosen so that the power series expansion of 
u[K/L/M](s) agrees with that of u(s) through the first N = K + L + M + 1 
terms. The constant terms in P and Q can be set arbitrarily to one without 
loss of generality since these choices do not alter the roots of the equation, so 
the total number of degrees of freedom isN = K + L + M + l. As true for 
ordinary Pade approximants, the coefficients of the polynomials can be com- 
puted by solving a matrix equation and the most accurate approximations are 
obtained by choosing the polynomials to be of equal degree, so-called "diago- 
nal" approximants. Because the power series begins with s 4 , the method was 
applied to u = u(x[s]) / (xqS 4 in the coordinate s; the eigenparameter is then 
estimated from 

x ss u(s = -1)~ 1/5 (18) 
The quadratic equation actually yields two approximations: 



± , -Q(s) ± x /Q(s)i-4P(s)R(s) 



u±[K/L/M](s) = - (19) 



as illustrated for K = L = M = 16 in Fig. [TJ the root obtained from applying 
the minus sign in front of the radical is spuriously negative on much of the 
interval s E [—1,0]. Substituting u = x^u^lK / L / M] into the differential 

equation yields the "residual function, p ± = u ss /xl — Jxq(s + l)u. Plotting 
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the residual functions as in Fig. [2] shows that the "physical" root is the plus 
sign in ( TT9|) ; the residual function is tiny for the plus root but 0(1) for the 
negative root. One warning is that, like ordinary Pade approximants, Hermite- 
Pade computations are rather ill-conditioned. We therefore employed Maple 
so that much of the work was done in exact rational arithmetic, and floating 
points portions were calculated using 60 decimal digits of precision; we did 
not investigate ill-conditioning. 

The first surpise is that both solution branches of the Hermite-Pade converge to 
almost the same values as s f> -1. We therefore list the errors of the diagonal 
Shafer approximants for xq from both solutions in Table [6j The table also gives 
the error of the power series up to order N from whence the Hermite-Pade 
approximations were obtained. We noticed that the errors of the two branches 
are roughly equal and opposite; we therefore also list the error in the average 
of the two roots of the quadratic. Each of the roots of the quadratic is a much 
better approximation than the power series. 

In other applications of Hermite-Pade approximations, the roots converge to 
separate modes, such as one root to the ground state eigenvalue and the other 
to the first excited state. Here, however, both roots converge to the unique 
eigenparameter and their average is extraordinarily accurate. We are unaware 
of another profitable use of averaging in Hermite-Pade approximations. 

The second surprise is that when K > 17, the diagonal approximant develops 
a spurious singularity in the middle of the spatial interval for both branches as 
illustrated in Fig. |3j Except in a narrow neighborhood of the singularity, the 
residual function for one branch is very tiny as illustrated in Fig. HJ It has long 
been known that ordinary Pade approximants may develop similar spurious 
singularities even at high order and even in spatial regions where the approxi- 
mants have started to converge; the difficulty is not due to roundoff error, but 
rather to a near-coincidence of zeros in the numerator and denominator of the 
rational function which is the Pade approximant [121H3"] . Similarly, the three 
coefficients of the Hermite-Pade quadratic and also the discriminant, which is 
the argument of the square root in the approximations, all have nearly coin- 
cident zeros as illustrated in Fig. The branches also switch identities at the 
singularities so that the physical solution is given by the "plus" branch on one 
side of the singularity and by the "minus" branch on the other side. 

This difficulty is not as serious as it seems. Both branches and their average 
continue to give superb approximations to the eigenparameter xq. The ordi- 
nary power series converges exponentially fast in the middle of the spatial 
interval where the Hermite-Pade fails. 

The third surprise is that when the order of the approximant is thirty or 
larger, both roots predict complex-valued x . When K = L = M = 30, for 
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example, the two branches of the quadratic give errors in xq of 0.199 x 10 -23 + 
i0.166 x 10~ 15 . The imaginary parts cancel when the average is taken, leaving 
the extraordinarily small and real valued average error of about 0.2 x 10 -23 ! 
Thus, in spite of the three surprises, which are also complications, and need 
for high precision arithmetic at high order, the Hermite-Pade acceleration of 
the right endpoint Taylor series is very successful. 



4-3 Power-series for a modified equation 



We can improve the results by removing the singularity through a conve- 
nient transformation of the differential equation (jSJ). Instead of the variables 
discussed avove it is more convenient for our aims to choose x = xqz 2 and 
v(z) = \Ju{xqz 2 ) so that the differential equation becomes 

zv(z)v"(z) + zv'(z) 2 - v(z)v'{z) - 2x 5 Q /2 z A v(z) = (20) 



In this way the domain size xq appears as a sort of eigenparameter xq. Fig [6] 
compares the two functions u(x) and v (z). 

We can thus obtain xq and u' from the partial sums 



M 

v M {z)=x 5 /2 t 2 J2cj t: ' 

j=0 



5/2 2 /l 10 149 2 . 

= i 'T - + — t + 1 2 + . . . , t = z - 1 21 

V3 21 588 J ' V ' 

and the equations fju(0) = 1 and u' Q = v'^(Q)/xq. Results can be considerably 
improved by means of Pade approximants [J,K](t). Here we choose diagonal 
[M/2, M/2] and near-diagonal [(M — 1)/2, (M + l)/2] ones because experience 
and theory suggest that they are usually the most accurate [4"2H4"1] . Tables [7] 
and [8] show the parameters u' and xq calculated by means of the partial sums 
(|2ip and their Pade approximants. The rate of convergence is remarkable for 
a power-series approach to a nonlinear problem. 



5 Chebyshev Pseudospectral Method 



We have also written a program that solves the Thomas-Fermi equation using 
the Chebyshev pseudospectral method and Newton iteration. Because it is so 
similar to an earlier for the Lane-Emden equation, we omit the details [17]. 
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5. 1 Calculations of large order 



The solution is approximated in the form 



oo 



v ( z ) = a nT n (2z ~ 1) 



(22) 



n=0 



The Chebyshev polynomials can be conveniently evaluated by T n (z) = cos(n arccos(z)) 
or by the usual three term recurrence relation [16] . 

Newton's iteration, which was used to solve the system of quadratic equations 
for the unknowns (x , a , ai, . . . a at), requires and initialization or "first guess". 
We found that a first guess for x £ [2, 6] combined with the lowest term in 
the power series, v ~ x^ 2 (z — l) 2 /3, suffice to give rapid convergence without 
underrelaxat ion . 

The Chebyshev coefficients converge geometrically as illustrated in Fig. [7] with 
a n oc exp(— nfi) where \i is between 1.4 and 1.42. All the Chebyshev coefficients 
larger than 10 -15 listed in Table [9J The Chebyshev series converges much more 
rapidly than the power series about the right endpoint as illustrated in Fig. |HJ 

Based on the trends in the Chebyshev calculations, we believe the following 
results are accurate to all 50 decimal places shown. 

x = 3.06885718281479942624073100623167158584582595057745 (23) 

u ' = -0.93896688764395889305505340187460180383289370739437 (24) 

5.2 Semi- Analytical Solutions: Chebyshev Polynomial Methods for Small N 

The most efficient spectral representation is one that incorporates all three 
boundary conditions into the approximation in a form independent of the 
spectral coefficients d n : 



v = (z- l) 2 {1 + d 1 [T 1 (2z - 1) + T ] + d 2 [T 2 (2z - 1) + T x (2z - 1) + . . .} (25) 



where we have used the identity sign(T n (— 1)) = — sign(T n+1 (— 1)) for all n. 
It is convenient, to eliminate fractional powers in the unknowns, to define the 
new unknown 




(26) 
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The pseudospectral method with collocation at the interior points Xj = (1 + 
cos(irj/(N + 2)), j = 1, . . . (N + 1) yields a system of coupled polynomial 
equations, quadratic in the unknowns (A, d x , d 2 , . . . d^). 

For N = 1, for example, substitution of v — (z — l) 2 (1 + 2d x z) into the 
differential equation yields the residual function 



resid(z; A, di) = z v v" + z(v ') 2 — v v' — A z 4 v 

= 2 - 2 A z 7 a + 4 A z 6 a - 2 A z 5 a - 6 z 2 - 64 z 3 a + 36 a z 2 - 120 z A a 2 Q 
+ 96 2 3 a 2 - 24 a 2 ^ 2 + 30 z A a + 48 z 5 a 2 - A z 6 + 2 A z 5 - A / + 4 - 2 a 

(27) 

The coupled system of two equations in two unknowns is the pair of equations 
rem?(l/4; A, d\) = and resz(i(3/4; A, c?i) = 0: 



27 81 j 27 j 2 9 . 9 . , n 

di d\ A A d\ = 

16 128 64 4096 8192 

5 95 , 27 , 2 81 . 243 . , n 

1 tZH di A A = (28) 

16 128 64 4096 8192 v ; 

The resultant of the two equations with elimination of d\ is 

2705319 1643895 , 51182361 l2 1240029 ^ 

4194304 268435456 68719476736 2199023255552 1 ; 



This, with substitution of each A root in turn back into one of the residual 
conditions to determine d ly yields the three solutions: 

{A = 34.33616, d x = 1.31544}; {A = 1311.8, d x = -.66653}; {A = -25.3931, d x = -2.8724}(30) 



The first solution is graphed in figure [HI the other two are spurious. 

For small N, the Maple "solve" command will find all the finite solutions 
to the system, here 2 N+1 — 1 in number. In this work, the physical solution 
was identified as that solution with the eigenparameter A closest to known 
values found by other means. When no such a priori information is available, a 
tedious but reliable procedure is to substitute each solution into the differential 
equation to calculate the residual, and accept the solution with the smallest 
residual norm. If multiple solutions are suspected, a good strategy is to use 
each small- Absolution as a first guess for the Newton pseudospectral code at 
higher- resolut ion . 

Maple's system solver is very slow and failed for A^ = 5; in contrast, the 
Newton/collocation program needed only 60 seconds to calculate for A^ = 100 
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in 100 decimal place arithmetic. However, the Maple code that exploits the 
" solve" command is much shorter and is given in its entirety in Table [TO] 

The numerical results, including coefficients up to and including d^, are in 
Table [11] and Fig. [9j Note that only the N — 1 approximation, obtained by 
solving a pair of equations for (A, di), is graphically distinguishable from the 
exact v(z). The success of such low truncations emphasizes the smoothness of 
the Thomas-Fermi solution in the transformed coordinate z and unknown u. 
it also illustrates that spectral methods can often be used in a semi-analytical 
mode as described further in [15] and Chapter 20 of [36]. 



6 Conclusions 

We have shown that the PHM converges for two Thomas- Fermi equations that 
exhibit quite different boundary conditions. In the case of neutral isolated 
atoms we have u(x — > oo) = and u'(x — > oo) = while, on the other 
hand, u(x ) = and u'(x ) = apply to the neutral atoms in very strong 
magnetic fields. If the success of the PHM depends on the existence of a 
movable singularity that the approach can push to infinity in the form of a 
zero of the denominator of the Pade approximant |23) . then it seems that 
two different kinds of singularities are involved in the problems just outlined. 
According to Abbasbandy and Bervillier [23] , in the former case such singular 
point is precisely the pole x s (see Eq. However, in the latter case the only 
candidate appears to be the zero/branch-point x q (see Eq. (TTUT) ) . Besides, in 
this case the rate of convergence of the PHM is remarkably larger. 

Present PHM calculations of the critical slope at the origin are quite accurate 
but the generation of analytic Hankel determinants of dimension as large as 
D = 60 is time consuming. However, the rate of convergence is commonly so 
great that one can obtain reasonably accurate results from determinants of 
relatively small dimension. This feature of the PHM was already exploited by 
Abbasbandy and Bervillier [23] to estimate the parameters of the conformal 
mapping used in the analytical continuation of the power series. 

We have also shown that nowadays available computer algebra software like 
Mahematica enable us to obtain the solution to a nonlinear differential equa- 
tion quite efficiently and accurately by means of suitable built in commands. 
Such results are shown in tables HI H] and [5j However, we found the PHM more 
convenient for the accurate calculation of the slope at the origin as shown in 
tables □ and H 

In the case of the Thomas- Fermi equation for a neutral atom in a very strong 
magnetic field we have shown that the power-series expansion about the zero 
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of multiplicity two at x is a suitable way of obtaining reasonably accurate 
approximants to the solution over the entire physical interval < x < Xq. 
The results can be slightly improved by means of Pade approximants and 
considerably improved by means of Hermite-Pade approximations, even when 
the latter exhibit a singular point or yield complex results at large orders 
of approximation. If we remove the singularity at origin by means of a suit- 
able transformation both the power series and its Pade approximants lead to 
remarkably more accurate results. 

Finally, we have proved once more that the Chebyshev polynomials are by far 
the most accurate and efficient way of solving this type of equations. In this 
application we resorted to another computer-algebra software, Maple. 
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Ref. J22] 
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Ref. [23] 
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Ref. El 
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Hyper asymptotics, Acta Applicandae Mathematicae 56 (1999) 1-98. 
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Table 2 

Critical slope at origin for the Thomas-Fermi equation for neutral atoms in strong 
magnetic fields 



Source u' 

Ref. [8] -0.93896594 

Ref. p] -0.938966887644 

Ref. [22] -0.93896688764395889306 

Present (PHM) -0.9389668876439588930550534018746018038328937073944 
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Table 3 

Solution to the Thomas-Fermi equation for neutral atoms 

x u(x) 

1. 

10 0.024314292988303278642 

20 0.0057849411914221299709 

30 0.0022558366161877068501 

40 0.0011136356388322203339 

50 0.0006322547829808588775 

60 0.00039391136668444560951 

70 0.00026226529980867263992 

80 0.00018354575974054805106 

90 0.00013354582895236378464 

100 0.00010024256813825348287 

110 0.000077192183912095346003 

120 0.000060724454047393954368 

130 0.000048642170610259129166 

140 0.000039574139447381720342 

150 0.00003263396444616844078 

160 0.000027231036946022553785 

170 0.000022961351006569413655 

180 0.000019542101672232006661 

190 0.000016771248040459921349 

200 0.000014501803496663821593 

210 0.000012625078687061407726 

220 0.00001105951450017485918 

230 9.743090050747086742 x 10~ 6 

240 8.628066979117904004 x 10~ 6 

250 7.6772907644292389823 x 10~ 6 

260 6.8615483805952773187 x 10~ 6 

270 6.1576544141466258996 x 10~ 6 

280 5.5470471164203299498 x 10~ 6 

290 5.0147463892445580293 x 10~ 6 

300 4.548571953477564028 x 10~ 6 

310 4.13855079r§736951774 x 10~ 6 

320 3.7764638005770180322 x 10~ 6 



Table 3 
Continuation 



X 




u(x) 






360 


2. 


.6868954789918090317 x 10" 


-6 


370 


2. 


.4819436135216654243 x 10" 


-6 


380 


2, 


,2973543392649800492 x 10" 


-6 


390 


2, 


.1306570419316347382 x 10" 


-6 


400 


1, 


.979732628107016001 x 10~ 6 




410 


1, 


.8427564849149959806 x 10" 


-6 


420 


1 


.7181517839000305473 x 10" 


-6 


430 


1 


.6045510643415443124 


x 10" 


•6 


440 


1, 


.5007644808004633506 


x 10" 


-6 


450 


1. 


.4057534397292700168 x 10" 


-6 


460 


1, 


.3186086183010787724 x 10" 


-6 


470 


1. 


.2385315617209843353 x 10" 


-6 


480 


1 


.164819216583818838 x 10~ 6 




490 


1, 


.0968508828337335859 x 10" 


-6 


500 


1. 


.0340771681908360586 x 10" 


-6 


510 


9. 


,7601060361700760275 


x 10" 


-7 


520 


9. 


.2221764588916452024 


x 10" 


-7 


530 


8. 


.7231183937133146868 x 10" 


-7 


540 


8. 


.259479517181917443 x 10~ 7 




550 


7. 


.8281693038804873437 


x 10" 


-7 


560 


7, 


.4264155192802616799 


x 10" 


-7 


570 


7. 


.0517266033530000494 


x 10" 


-7 


580 


6. 


.7018590438324510541 


x 10" 


-7 


590 


6. 


.3747890207491997214 


x 10" 


-7 


600 


6. 


.0686876966498027289 


x 10" 


-7 


610 


5. 


.7818996334617716871 


x 10" 


-7 


620 


5. 


.5129238989702892144 


x 10" 


-7 


630 


5, 


.2603974914209643451 


x 10" 


-7 


640 


5. 


.0230807665120921283 


x 10" 


-7 


650 


4. 


.7998445982454301833 


x 10" 


-7 


660 


4, 


.5896590452695361208 


x 10" 


-7 


670 


4. 


.39158332§@887189732 


x 10" 


-7 


680 


4. 


.2047569472008565958 


x 10" 


-7 



Table 3 
Continuation 



X 


u{x) 




710 


3.7042134434543470637 x 10" 


-7 


720 


3.5551266392993614555 x 10" 


-7 


730 


3.4139433122891557142 x 10" 


-7 


740 


3.2801461578314272486 x 10" 


-7 


750 


3.1532580023535960183 x 10" 


-7 


760 


3.0328382170354037711 x 10" 


-7 


770 


2.9184794904453079974 x 10" 


-7 


780 


2.809804924949050881 x 10~ 7 




790 


2.7064654197531601863 x 10" 


-7 


800 


2.6081373047370408811 x 10" 


-7 


810 


2.5145202065854619574 x 10" 


-7 


820 


2.4253351127307245963 x 10" 


-7 


830 


2.3403226197770901371 x 10" 


-7 


840 


2.2592413435921192214 x 10" 


-7 


850 


2.1818664751032822416 x 10" 


-7 


860 


2.1079884666443520651 x 10" 


-7 


870 


2.0374118362107619892 x 10" 


-7 


880 


1.9699540776667432384 x 10" 


-7 


890 


1.9054446664164008206 x 10" 


-7 


900 


1.8437241512360732429 x 10" 


-7 


910 


1.7846433238508584078 x 10" 


-7 


920 


1.7280624596140222668 x 10" 


-7 


930 


1.6738506202412198291 x 10" 


-7 


940 


1.621885016469136105 x 10" 7 




950 


1.5720504224086543073 x 10" 


-7 


960 


1.5242386363458728059 x 10" 


-7 


970 


1.4783479864166607799 x 10" 


-7 


980 


1.4342828759880769669 x 10" 


-7 


990 


1.3919533628465534453 x 10" 


-7 


1000 


1.3512747734085860175 x 10" 


-7 
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Table 4 

First derivative of the Solution to the Thomas-Fermi equation for neutral atoms 



X 




v!{x) 





-1 


.5880710226113753127 


10 


-0. 


.0046028819145712279321 


20 


-0. 


.00064725433367546009732 


30 


-0. 


.00018067000619215999232 


40 


-0. 


.000069668028734127283349 


50 


-0. 


.000032498901894147568391 


60 


-0. 


.000017197699996624380177 


70 


-9 


.9565334520339929033 x 10 - ' 


80 


-6 


.1661955042322463779 x 10"' 


90 


-4. 


.0244737401132216063 x 10 - ' 


100 


-2 


.7393510441626592228 x 10"' 


110 


-1. 


,9299022618702733326 x 10 - ' 


120 


-1. 


,3992583313630533608 x lO - ' 


130 


-1. 


.0395157141360444629 x 10 - ' 


140 


-7. 


.8856446105551487333 x 10" 


150 


-6 


.0913995255007820982 x 10~ 


160 


-4. 


.7807415119622515851 x 10 - 


170 


-3 


,8051134233773750717 x 10" 


180 


-3 


,0666432505550168717 x 10 _ 


190 


-2 


.4992901359390774381 x 10" 


200 


-2. 


.0575322915697839603 x 10" 


210 


-1. 


.7093868083213874424 x 10" 


220 


-1. 


.4319926122797382022 x 10" 


230 


-1. 


.2087522194709907058 x 10 _ 


240 


-1. 


,0274430678617277247 x 10" 


250 


-8. 


.7894678848684294245 x 10- 


260 


-7 


.5637915959323113671 x 10- 


270 


-6 


.5448528599341097551 x 10~ 


280 


-5. 


.6921313039625962206 x 10~ 


290 


-4. 


.9740860613513071015 x 10~ 


300 


-4. 


,3659496739887766774 x 10 _ 


310 


-3 


.84811438^)621989067 x 10" 


320 


-3 


.4049389150220980733 x 10~ 



Table 4 
Continuation 



X 




u'(x) 






360 


-2. 


.1605807533839970561 


x 10" 


-8 


370 


-1. 


,9432625416385162874 


x 10" 


-8 


380 


-1. 


.7526290456713803171 


x 10" 


-8 


390 


-1. 


.5848392596773086629 


x 10" 


-8 


400 


-1. 


.4366822938881119184 


x 10" 


-8 


410 


-1. 


.3054622274557951893 


x 10" 


-8 


420 


-1. 


.1889056354032339405 


x 10" 


-8 


430 


-1. 


.0850874762054582989 


x 10" 


-8 


440 


-9 


.9237153282410934663 


x 10" 


-9 


450 


-9 


.0936175646022167788 


x 10' 


-9 


460 


-8. 


.3486284112277083368 


x 10" 


-9 


470 


-7. 


.6784786402098636917 


x 10" 


-9 


480 


-7. 


.07431704734312557 x 


io- 9 




490 


-6 


.5284907212008860465 


x 10" 


-9 


500 


-6 


,0343633675250736044 


x 10" 


-9 


510 


-5. 


.5861638814315092358 


x 10' 


-9 


520 


-5. 


.178858949068817509 > 


< 10- 9 


530 


-4. 


.8080478711684610573 


x 10" 


-9 


540 


-4. 


,4698732556423909026 


x 10" 


-9 


550 


-4 


,1609451793597615307 


x 10" 


-9 


560 


-3 


,8782777480051922656 


x 10" 


-9 


570 


-3 


.6192350163757433335 


x 10" 


-9 


580 


-3 


.3814850215456566458 


x 10' 


-9 


590 


-3 


.1629599078554295417 


x 10" 


-9 


600 


-2 


.9618225460409763933 


x 10" 


-9 


610 


-2. 


,7764375594882403909 


x 10" 


-9 


620 


-2. 


,6053463478912139361 


x 10" 


-9 


630 


-2. 


.4472456671467371541 


x 10" 


-9 


640 


-2. 


.3009688041347206372 


x 10" 


-9 


650 


-2. 


.1654695131200772501 


x 10" 


-9 


660 


-2 


,0398080382166108813 


x 10" 


-9 


670 


-1. 


,92313912^)028323146 


x 10" 


-9 


680 


-1. 


.8147008592751772428 


x 10" 


-9 



Table 4 
Continuation 



X 




u'(x) 






710 


-1 


.5322179719643453493 x 10 


-9 


720 


-1 


.4504531300856574909 x 10 


-9 


730 


-1 


.3740749690000988024 x 10 


-9 


740 


-1 


.3026628314027827248 x 10 


-9 


750 


-1 


.2358340874513047305 x 10 


-9 


760 


-1. 


.1732399946130836533 x 10 


-9 


770 


-1 


.1145622752274576869 x 10 


-9 


780 


-1 


,059510454133309805 > 


< 10~ : 


3 


790 


-1 


.0078187598866428379 x 10 


-9 


800 


-9. 


.592438634509292608 x 10 - 


10 


810 


-9. 


.135628191336062491 x 10 - 


10 


820 


-8. 


.7057109477720637091 


x 10 


-10 


830 


-8. 


.3008081231100632109 


x 10 


-10 


840 


-7. 


.9191919524473132439 


x 10" 


-10 


850 


-7. 


.5592725659439472365 


x 10 


-10 


860 


-7. 


.2195855931622667064 


x 10 


-10 


870 


-6. 


.8987806912593494372 


x 10" 


-10 


880 


-6. 


.5956115226598912448 


x 10 


-10 


890 


-6. 


,3089267086832537466 


x 10 


-10 


900 


-6. 


.0376616555034554537 


x 10" 


-10 


910 


-5. 


.7808313244166361146 


x 10" 


-10 


920 


-5. 


.5375229966408037985 


x 10" 


-10 


930 


-5. 


.306890430165250796 x 10" 


10 


940 


-5. 


.0881496461401375725 x 10 


-10 


950 


-4. 


,8805717870731175806 x 10 


-10 


960 


-4. 


.6834792695675137503 x 10" 


-10 


970 


-4. 


.4962431433435273583 x 10 


-10 


980 


-4. 


.3182774111530835327 x 10" 


-10 


990 


-4. 


.1490358684843996374 x 10 


-10 


1000 


-3 


.9880107279139957822 x 10 


-10 
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Table 5 

Solution to the Thomas-Fermi equation for neutral atoms in strong magnetic fields 

x u(x) 

1 

0.1 0.9069294578 

0.2 0.8167818236 

0.3 0.7306447995 

0.4 0.6491710501 

0.5 0.5727799968 

0.6 0.5017304847 

0.7 0.4361596667 

0.8 0.376107784 

0.9 0.3215357471 

1 0.2723385646 

1.1 0.2283561701 

1.2 0.1893825135 

1.3 0.1551734414 

1.4 0.1254536975 

1.5 0.09992326627 

1.6 0.07826321273 

1.7 0.06014112641 

1.8 0.04521625071 

1.9 0.03314435629 

2 0.02358240435 

2.1 0.0161930348 

2.2 0.01064890732 

2.3 0.006636917141 

2.4 0.003862303605 

2.5 0.002052666026 

2.6 0.0009618988483 

2.7 0.0003740561223 

2.8 0.0001071536822 

2.9 0.00001691610305 

3 0.0000004744716194 
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Table 5 
Continuation 



x u(x) 

3.1 0.00000002016132394 

3.2 0.000006420928387 

3.3 0.00006280071944 

3.4 0.0002680829463 

3.5 0.0007805068545 

3.6 0.001821117342 

3.7 0.003677230211 

3.8 0.006705874884 

3.9 0.01133721636 

4 0.01807795804 

4.1 0.02751472686 

4.2 0.04031744196 

4.3 0.05724266818 

4.4 0.07913695529 

4.5 0.106940164 

4.6 0.1416887798 

4.7 0.1845192145 

4.8 0.236671098 

4.9 0.2994905589 

5 0.3744334963 

5.1 0.4630688413 

5.2 0.5670818119 

5.3 0.6882771577 

5.4 0.8285823989 

5.5 0.9900510563 

5.6 1.174865876 

5.7 1.385342047 

5.8 1.623930412 

5.9 1.893220673 

6 2.195944591 
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Table 6 

Errors in Eigenlength xq for the magnetic case, obtained from right endpoint power 
series and Hermite-Pade approximation for u{x) 



K (= L = M) 


N 


Taylor series 


first root 


second root 


average of roots 


2 


7 


0.0039 


0.0023 


0.044 


0.023 


5 


16 


0.00076 


0.00093 


0.000052 


0.00043 


10 


31 


0.174 x 10" 3 


0.265 x 10~ 6 


0.32 x 10~ 6 


0.31 x 10~ 7 


15 


46 


0.69 x 10~ 4 


0.11 x 10~ 8 


0.11 x 10~ 8 


0.32 x lO" 11 


20 


61 


0.35 x 10~ 4 


0.15 x 10~ n 


0.15 x lO" 11 


0.45 x 10~ 15 


25 


76 


0.21 x 10~ 4 


0.11 x 10~ 12 


0.11 x 10~ 12 


0.39 x 10~ 18 



Table 7 



v-power series approximation of initial slope u' Q ; correct digits in boldface 



M 


Method 


u' 


Number of correct digits 


10 


power series 


-0.9364 


2 




Pade 


-0.9327 


2 


20 


power series 


-0.93879 


3 




Pade 


-0.9389645 


5 


40 


power series 


-0.938966834 


7 




Pade 


-0.93896688760 


9 


60 


power series 


-0.938966887635 


10 




Pade 


-0.9389668876439554 


14 


80 


power series 


-0.9389668876439574 


14 




Pade 


-0.9389668876439588927 


17 




PHM [22] 


-0.93896688764395889306 


20 
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Table 8 



Eigenparameter xq approximations from the -u-power series 



M 


Method 


Xq 


10 


power series 


3.068806 




Pade 


3.06877 


20 


power series 


3.0688560 




Pade 


3.068857176 


40 


power series 


3.06885718271 




Pade 


3.068857182814799451 


60 


power series 


3.0688571828147917 




Pade 


3.06885718281479942624073139 


80 


power series 


3.0688571828147994255 




Pade 


3.068857182814799426240731006231672626130 
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Table 9 

Chebyshev coefficients of the solution to the magnetic case of the Thomas-Fermi 



equation 



degree n 


a n 





5.48135788388634344 x 10- 


01 


1 


-5.58881613523477153 x 10- 


-01 


2 


-6.42705283349619230 x 10- 


-02 


3 


5.70880541499120995 x 10- 


02 


4 


1.61861107728602653 x 10- 


02 


5 


1.78344934649533460 x 10- 


03 


6 


-4.97638953814069498 x 10- 


-05 


7 


9.88568187276911699 x 10- 


06 


8 


-1.57314041941070364 x 10- 


-06 


9 


2.17031601079801803 x 10- 


07 


10 


-3.17277060187282838 x 10- 


-08 


11 


6.73733273124207577 x 10- 


09 


12 


-1.91849909540107163 x 10- 


-09 


13 


5.43530178894964329 x 10- 


10 


14 


-1.37539658139417754 x 10- 


-10 


15 


3.12411846416475022 x 10- 


11 


16 


-6.65185565734159168 x 10- 


-12 


17 


1.41167490452475800 x 10- 


12 


18 


-3.15485896233555226 x 10- 


-13 


19 


7.52433060541288362 x 10- 


14 


20 


-1.85848288796828985 x 10- 


-14 


21 


4.57976772880843444 x 10- 


15 


22 


-1.10484891112084068 x 10- 


-15 
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Table 10 

Maple Code to Compute a Chebyshev Pseudospectral Approximation for Small N 
restart; with(orthopoly); N:=4; # Define approximation in next 2 lines; 

vv:=l; for n from 1 to N do vv:=vv + d[n]*(T(n,2*x-l)+T(n-l,2*x-l)); od: 

v:= (x-l)*(x-l)*w; # derivatives of v; vx:=diff(v,x); vxx:=diff(v,x,x); 

resid:= x*v*vxx + x*vx*vx -v*vx - lambda* x**4 *v; # ODE residual; 

ta[0]:=evalf( Pi*(0+l)/(N+2) ); xa[0]:=evalf((cos(ta[0])+l)/2); # Chebyshev grid points "xa"; 
resida[0]:= evalf(subs(x=xa[0],resid)); 

varset:= lambda; eqset:= resida[0]; # initialize set of unknowns & the set of pointwise residuals; 
for j from 1 to N do ta[j]:= evalf( Pi*(j+l)/(N+2) ); 

xa[j]:= evalf( (cos(ta[j] ) + l)/2 ); resida[j]:= evalf(subs(x=xa[j],resid)); # array of ODE residuals ; 

varset:= varset union d[j]; eqset:= eqset union residafj]; # update variable & equation sets; od: 

solutionvector:=solve(eqset,varset); assign(solutionvector[l]); # solve the polynomial system; 

# {resida[j]=0, j=0, N} in unknowns {A, di, . . . d^} ; 

xi:= evalf( (lambda/2)**(2/5) ); sigma:= evalf( subs(x=0,diff(v,x,x)/xi) ); 



Table 11 

Results and Errors from Low Order Collocation 



Name 


N = 1 


N = 2 


N = 3 


N = 4 


Exact 




3.11809 


3.07417 


3.0687055 


3.068876456 


3.0688571 


xo error 


0.049 


0.0053 


-0.00015 


0.19274 x 10-4 




u' 


-2.73 


-.281 


-.89047 


-.9237607 


-0.9389669 


u' error 


-1.79 


0.66 


0.048 


0.015 




di 


1.315 


1.8211 


1.8837159 


1.882043824 


1.8825071389 


d 2 




.29256 


.34172159 


.3402298497 


.34061038 


d 3 






0.027764 


0.0266875 


.0269894965 


C?4 








-0.000662 


-.0004613479537 
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ODE residuals K=L=M=16 Hermite-Pade 
s 

-0.8 -0.6 _ 94——^ _ -n? 




Thomas-Fermi magnetic 



,1e2 
1. 

e-1 

e-2 

e-3 

e-05 

e-06 

e-07 

e-08 

e-09 

e-10 

e-1 1 

e-1 2 

e-1 3 

e-1 4 

e-1 5 

e-1 6 

e-1 7 

e-18 

e-1 9 

e-20 

e-21 

e-22 

e-23 



Fig. 2. The K = L = M = 16 Hermite-Pade ODE residuals, that is, the result of 
substituting the approximation into the Thomas-Fermi differential equation. 
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K=L=M=20 Hermite-Pade: both branches 




Fig. 3. The K = L = M = 20 Hermite-Pade approximations to u in the coordinate 
s = z/xq — 1. 



ODE residuals of Hermite-Pade for Thomas-Fermi magnetic 
K=L=M=20 




Physical branch 



Fig. 4. The K 

S = z/xq — 1. 



L = M = 20 Hermite-Pade ODE residuals in the coordinate 
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Hermite-Pade: P, Q, R & Discriminant K=L=M=20 



.le-1-. 




-0.601 -0.600 -0.599 -0.598 -0.597 -0.596 



Fig. 5. The K = L = M = 20 Hermite-Pade approximant: Plots of the coefficient 
polynomials in the quadratics equation P(f[K/L/M}) 2 + Q f[K/L/M] + R = 
and also the "discriminant" , D = Q 2 — APR, which is the argument of the square 
root in the solutions to the quadratic. 




Fig. 6. Left: the solution to the Thomas-Fermi problem or strong magnetic field. 
Right: the transformed unknown v(z) in the rescaled coordinate space z. 
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v,v1, Loo error=0.063904 



v,v2, Loo error=0.0053463 





v, v3, Loo error=0.0001 231 9 v, v4, Loo error=3.5303e-005 

1 r— =^ 1 1 





Fig. 9. Comparisons of v(z) [solid black] with N = 1 to N = 4 Chebsyhev pseu- 
dospectral approximations [dashed red]. The green dotted lines, barely visible at 
the bottom, are the errors. 
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